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A Two  Dimensional  Global  Climatic  Model 


1 . Summary 

The  purpose  of  this  project  is  to  develop  a simple  global 
climatic  model  which  may  be  used  to  study  the  distribution  of  climaces 
on  the  earth  and  their  causes,  the  factors  that  influence  climatic 
change,  and  the  possible  influence  that  man  might  have  on  climate. 

The  model  structure  was  described  in  the  previous  Technical  Report 
for  this  project.  Briefly,  it  utilizes  a global  box  grid,  with  each 
box  having  dimensions  of  10°  of  latitude  by  10°  of  longitude.  The 
rates  at  which  humidity  in  the  atmosphere  and  temperature  and  the 
north-south  wind  or  current  speed  in  the  atmosphere  and  oceans  vary 
vertically  are  specified.  The  vertically-averaged  thermodynamic  energy 
equation,  the  equation  of  motion  for  the  surface  boundary  layer,  and 
the  hydrostatic  equation  are  then  solved  jointly  for  each  box  at  one- 
month  time  steps  to  determine  the  global  distribution  of  sea-level 
pressure,  temperature,  and  wind  speed  and  direction. 

A few  small  changes  have  been  made  in  the  parameterization  of  the 
model.  These  involve,  principally,  the  sea-level  pressure  field,  the 
divergence  of  the  zonal  transport  of  heat,  water  vapor,  and  potential 
energy  in  the  atmosphere,  and  the  transport  of  heat  by  ocean  currents. 

Each  change  is  discussed  briefly  in  the  following  section. 

Also  described  in  the  following  section  is  a proposed  hyorologic 
cyde  to  be  used  with  the  model.  As  it  stands  now,  the  model  already 
yields  values  for  the  storage  rate  of  water  vapor  in  the  atmosphere 
and  the  divergence  of  the  water  vapor  flux  by  the  atmospheric  circulation. 
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These  may  be  combined  to  give  the  difference  between  monthly  evaporation 
and  precipitation  rates  for  each  grid  box.  Over  the  oceans  the  evapor- 
ation rate  may  be  estimated  from  established  formulae  relating  it  to 
the  near— surface  wind  speed,  temperature,  and  relative  humidity. 

Over  land  evaporation  rates  depend  mainly  on  the  amount  of  water  in 
the  soil  and  the  potential  evapotranspiration  rate  and  may  be  determined 
from  the  water  balance  equation  for  a soil  column  using  a simple  book- 
keeping approach.  As  by-products  here,  the  soil  moisture  content  and 
runoff  rates  are  obtained. 

Knowing  the  evaporation  rate,  the  monthly  precipitation  may  be 
obtained  directly  (with  the  necessary  constraint  that  it  may  never  be 
less  than  zero).  Depending  on  the  surface  temperature,  a prescribed 
fraction  of  this  precipitation  will  fall  as  snow,  adding  to  that  which 
may  be  already  present  on  the  ground  or  ice.  Sea  ice  forms  if  the  air 
temperature  is  at  or  below  freezing  and  decreasing.  It  melts  from 
below  when  the  temperature  is  rising.  Both  snow  and  ice  may  melt  at 
the  surface  when  the  temperature  is  above  freezing.  Part  of  the  snow 
melt  over  land  can  infiltrate  into  the  soil,  increasing  the  soil  moisture 
content.  The  fraction  of  each  box  covered  by  snow  and/or  ice  is  a simple 
function  of  their  respective  depths. 

The  relative  humidity  over  water,  snow,  or  ice  may  be  expressed 
as  the  mid-ocean  value  (prescribed  as  a function  of  latitude)  modified 
by  the  zonal  advection  of  wetter  or  drier  air  from  upwind  regions. 

Over  land  the  surface  relative  humidity  depends  mainly  on  the  soil 
moisture  content.  To  complete  the  cycle,  cloud  cover  can  be  empirically 
estimated  from  the  precipitation  and  evaporation  rates  and  the  surface 
relative  humidity. 
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With  the  incorporation  of  a hydroloyic  cycle,  the  only  prescribed 
variable  in  the  model  (other  than  the  mid-ocean  relative  humidity) 
is  the  extraterrestrial  solar  radiation.  Everything  else  can  vary. 

With  such  a system,  incorporating  numerous  feedbacks,  both  positive  and 
negative,  it  is  possible  that,  without  proper  coupling,  one  process  may 
dominate  and,  after  a sufficiently  long  run  of  the  model,  produce 
nonsense  results.  Therefore,  it  will  be  necessary  to  make  one  or  more 
runs  of  about  100  model  years  to  check  the  stability  of  the  model. 

Any  further  modifications  or  applications  will  depend  largely  on  the 
results  of  this  experiment.  Assuming  that  it  is  successful,  there  are 
several  useful  projects  that  can  be  carried  out.  These  include  comparing 
results  of  a particular  experiment,  for  example  that  involving  a change 
in  the  amount  of  carbon  dioxide  in  the  atmosphere,  with  those  obtained 
using  one  of  the  three-dimensional  general  circulation  models.  Tests 
on  the  sensitivity  of  the  model  to  changes  in  the  prescribed  initial 
conditions  and  studies  of  the  stability  of  the  Arctic  and  Antarctic 
continental  and  sea  ice  would  also  be  useful. 

2.  Model  revisions 

a.  sea-level  pressure 

Instead  of  specifying  the  meridional  pressure  gradient  averaged  over 
all  longitudes,  it  is  now  expressed  as  a function  of  the  corresponding 

sea-level  temperature  gradient.  This  is  similar  to  the  procedure  used 

in  Sellers  (1973).  However,  the  constant  of  proportionality  is  now 

assumed  to  be  inversely  proportional  to  the  zonally-averaged  surface 
friction  coefficient,  which  equals  4 x 10”2 * * * * 7  cm-1  over  land  and  i:e 
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and  2 > 10" 8 cm"1  over  water.  This  permits  asymmetr^  in  the  wind 
field  between  the  two  hemispheres,  with  much  stronger  surface  westerlies 
in  the  southern  hemisphere. 

b.  divergence  of  the  zonal  transport 

Because  the  model  contains  no  stratosphere,  the  zonal  wind 

increases  without  limit  to  the  top  of  the  atmosphere.  This  yields 

unrealistically  large  values  for  the  divergences  of  the  zonal  transport 

of  water  vapor,  heat,  and  potential  energy.  To  get  around  this,  the 

magnitude  of  each  divergence  is  constrained  and  not  permitted  to 

exceed  about  10  kly  mo-1 . This  is  done  by  not  allowing  the  magnitudes 

of  a „u  . a Uu  , and  a „u  (see  Sellers,  1973)  to  exceed  20,  1.4A, 

and  4.9A  cm  sec"1,  respectively,  where  A is  the  box  area  in  units  of 

lO^cm2  and  u is  the  surface  zonal  wind  speed.  The  restriction  on 
o 

the  magnitude  of  the  zonal  vapor  flux  divergence  should  be  kept  in 
mind,  since  it  will  probably  have  an  influence  on  the  precipitation 
field  obtained. 

c.  ocean-current  heat  transport 

In  determining  the  meridional  transport  of  heat  associated  with 
wind-driven  ocean  currents  what  is  needed  in  the  model  is  the  meridional 
component  vwq  of  the  component  of  the  surface  wind-driven  current 
(directed  45°  to  the  right  of  the  wind  V in  the  northern  hemisphere) 
in  the  direction  of  the  net  mass  transport  (45°  to  the  right  of  the 
surface  current  and  90°  to  the  right  of  the  wind) . Following  Adpm 
(1970,  pp.  415,  416),  this  implies  that 


(1) 


where  the  minus  sign  applies  to  the  northern  hemisphere  and  the  plus 
sign  to  the  southern  hemisphere. 


11 

o 

• 

o 
o 
o 
1— ' 
U1 

to 

to 

If  I 

r* 

when 

|v|  > 

600  cm  sec”1 

= 0.003728(1 

f V | 

when 

|V|  < 

600  cm  sec-1 

1 

Initial  runs  of  the  model  with  this  formulation  indicate  instabilities 
arising  near  the  equator  because  of  the  small  values  of  the  coriolis 
parameter  f there.  To  get  around  this,  it  may  be  necessary  to  let 
f = 10"4  sec”1  at  all  latitudes, 
d.  the  hydrologic  cycle 

The  hydrologic  cycle  subroutine  will  be  embedded  within  the 
main  program,  so  that  most  quantities  will  be  determined  more  than 
once  during  a given  month,  the  final  values  used  being  those  that 
give  satisfactory  convergence.  This  will  be  explained  more  clearly 
in  the  next  section.  At  any  rate,  if  it  seems  in  the  following 
discussion  that  some  quantities  are  being  used  that  have  not  been 
evaluated  yet,  this  is  really  not  the  case.  They  have  been,  but 
they  are  subject  to  change  before  the  calculations  are  completed. 

The  first  step  is  to  estimate  the  difference  between  the 
rates  of  evaporation  E and  precipitation  r from  estimated  or 
previously-determined  values  for  the  vapor  flux  divergence  and  the 
change  in  atmospheric  moisture  storage.  This  is  done  for  each  grid 
box.  Over  the  oceans  and  snow-covered  land  the  s'  .iace  relative 
humidity  RH  of  the  j-th  box  is  estimated  from 


RH  = RH'  = 0 . 5 ( RH  + e’/ej  / 


(2) 
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where  RH^  is  the  (specified)  mid-ocean  relative  humidity,  e' 
is  the  surface  vapor  pressure  of  the  upwind  box  and  e is  the 
saturation  vapor  pressure  of  the  j-th  box.  The  evaporation  rate  in 
cm  mo-1  may  then  be  determined  from 

E = E . = E'  = 0,004  I V I e (1  - RH)  . (3) 

pt  1 1 s 

Over  snow-free  land  at  temperatures  above  271K  the  potential 
evapotranspiration  rate  Ep  in  cm  mo-1  may  be  estimated  from 

E = 0.0016 (T  - 271) 3 , (4) 

p s 

where  Tg  is  the  surface  temperature.  This  relationship  was  derived 
from  data  given  by  Mcllroy  and  Angus  (1964)  for  Aspendale,  Australia. 
Actual  evapotranspiration  from  land,  which  may  be  partially  snow- 
covered,  is  given  by 

E = WiEp(l-As)/wk  + E'As  , (5) 

where  w^  is  the  soil  moisture  content  at  the  start  of  the  given 
month,  w^  is  the  critical  soil  moisture  content,  below  which  evaporation 
proceeds  at  less  than  the  potential  rate.  It  is  assumed  to  equal 
0.75  of  the  moisture  content  w at  field  capacity.  A is  the 
fraction  of  the  box  covered  by  snow.  This  is  obtained  from 

Ag  = 0.1  sd  (sd  < 10  cm)  (6) 

where  s^  is  the  existing  value  for  the  average  snow  depth  in  cm 

(water  equivalent) . In  (5)  when  w > w,  , the  ratio  of  the  two  is 

1 k 
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set  equal  to  1.0.  Knowing  the  evaporation  rate  from  (3)  and/or 
(5)  and  the  difference  E-r,  the  precipitation  rate  may  be  obtained 
directly,  with  the  condition  that  it  cannot  be  less  than  zero. 

In  order  to  estimate  cloud  cover  and  the  surface  albedo,  both 
needed  in  the  radiation  calculations,  it  is  now  necessary  to 
redetermine  the  snow  and  ice  coverage  and  the  soil  moisture  content. 
To  do  this,  it  may  first  be  assumed  that  the  snowfall  s equals  0 
when  T is  greater  than  275K  and  equals  r when  T is  less  than 
265K.  In  between 


s = (27.5  - 0.1  Tjr  . (7) 

It  is  further  assumed  that  the  potential  snow  melt  equals  7.5  E'. 

The  actual  snow  melt  s , of  course,  cannot  exceed  the  amount  of 

m 

snow  present  initially  plus  the  snowfall.  If  s^  is  the  average 
snow  depth  over  the  box  (or  over  ice  on  the  oceans)  at  the  beginning 
of  a given  month,  then  its  value  at  the  end  of  the  month  is 


s,  = s’  + s - s A over  land 
d d ms 


(8) 


and 

Sj  — si  + (s  - s_  (A  ) )A.  over  water  . (9) 

ad  m s i 

An  exception  to  (8)  is  made  in  Antarctica  and  Greenland.  The  depths 
of  these  two  ice  caps  are  held  constant. 

A^  in  (9)  is  the  fraction  of  the  box  covered  by  ice.  It  is 
related  to  the  existing  value  for  the  average  ice  thickness  id  in 
cm  by 


At  = 0.003  id  (id  < 333  cm) 


(10) 
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The  ice  thickness  can  change  by  melting  i^  or  freezing  i^  from  below 
or  by  melting  i^  from  above  in  snow-free  regions  at  a rate  equal  to 
or  less  than  the  potential  rate  7.5E1.  It  is  assumed  that 


and 


6(T  - T')  , 

s o 

6(t;  - Ts>  , 


i s=  0 when  T > T ' (11) 

f s s 

i =0  when  T < T'  < 275K  , (12) 
m s s 


where  T1  is  the  surface  temperature  at  the  start  of  the  given 
s 

month.  The  numerical  value  of  the  constant  was  selected  on  the 
basis  of  a statement  made  by  Wittman  and  Schule  (Fletcher,  1966, 
p.  219)  to  the  effect  that  one  season's  growth  of  ice  varies  from 
157  to  183  cm.  In  (11)  im  cannot  exceed  the  total  amount  of  ice 
that  can  be  melted.  The  ice  thickness  id  at  the  end  of  the  month 
is  given  by 


+ if  - (im  + ^ " smVAi 


(13) 


The  heat  required  to  melt  the  ice  from  below  and  that  released  by 

the  formation  of  snow  and  ice  should  be  included  in  the  energy 

conservation  equation.  The  energy  used  to  melt  surface  snow 

and  ice  is  already  included  in  the  evaporation  term. 

With  s,  and  i,  obtained  from  (8)  or  (9)  and  (13) , A and  A.  may 
d d ° 

be  redetermined  from  (6)  and  (10) . The  surface  albedo  ag  over 
water  may  then  be  estimated  from 


a = (0.8  A + 0.6(1  - A ) ) A . + (1  - A . ) a , 

S S S J.  a,  w 


(14) 


9 


where  is  the  albedo  of  open  water  (prescribed  as  a function  of 
latitude  and  time  of  year)  and  values  of  0.8  and  0.6,  respectively, 
are  used  for  the  albedos  of  snow  and  ice. 

Over  land  the  surface  albedo  depends  to  some  extent  on  the 
soil  moisture  content.  This  may  be  estimated  from  the  water 
balance  equation  for  a soil  column.  Following  Sellers  (1965) 
and  assuming  that  runoff  rQ  is  given  by 
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Having  determined  w and,  hence  w^,  which  becomes  w^  for  the  following 
month,  the  surface  albedo  for  land  may  be  obtained  from 


a = 0.8A  + (1  - A ) (0.3  - 0.25  w/w  ) . (21) 

s s s in 

Further,  the  surface  relative  humidity  and  total  potential  evapo- 
transpiration  are  given  by 

RK  = 0.9(1  - As)w/wm  + AsRH'  (22) 

and 

V = V1  • v + v (23) 

To  conclude  the  cycle,  it  is  assumed  that  cloud  cover  n is 
related  to  the  relative  humidity  and  the  precipitation  and  total 
potential  evapotranspiration  rates  by 


n = 


(0.75  + 0.25 


r - E 
r + E 


£^)rh  , 
pt 


(24) 


where  the  coefficients  have  been  determined  so  as  to  give  and 
average  cloud  cover  of  slightly  greater  than  0.5  under  present 
conditions.  In  the  model,  then,  cloud  cover  increases  with 
increasing  precipitation,  increasing  relative  humidity,  and 
decreasing  potential  evapotranspiration  (or  decreasing  temperature) . 


3.  Model  flow  diagram 

A simple  flow  diagram  will  be  used  to  show  how  the  hydrologic 
cycle  is  incorporated  into  the  model.  This  is  presented  in  Figure  1 
and  explained  in  the  following  paragraphs. 


FfciTK  *,  C-  - 
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The  first  step  is  to  estimate  or  specify  the  sea-level  tempera- 
tures for  each  grid  box  at  the  beginning  and  middle  of  a given 
month  (block  0} . The  calculations  are  started  with  January  and 
the  two  temperatures  are  assumed  equal.  The  sea-level  pressure, 
wind,  and  ocean  current  fields  corresponding  to  the  specified 
temperature  field  are  computed  next  for  the  whole  grid  (block  1) . 
Using  estimated  relative  humidities  (either  specified  or  those 
for  the  previous  month) , the  sea-lev '1  vapor  pressure  is  then 
determined  for  each  box  (block  2) . 

The  calculations  in  the  next  three  blocks  are  carried  out  at 
one  longitude  i at  a time.  First,  the  temperature  and  moisture 
dependent  parameters  are  evaluated  (block  3) . These  are  principally 
the  a coefficients  in  Sellers  (1973),  relating  vertical  integrals 
in  the  thermodynamic  energy  equation  to  their  sea-level  values. 

Next  (block  4)  the  components  of  the  hydrologic  cycle,  including 
cloud  cover  and  surface  albedo  are  computed.  This  requires  specified 
or  predetermined  estimates  for  ice  and  snow  thickness  and  the  soil 
moisture  content,  as  outlined  in  the  previous  section.  Then,  in 
block  5,  the  sea-level  temperature  in  each  10°  latitude  box  is 
determined  from  the  energy  equation. 

When  this  cycle  is  completed  for  all  longitudes,  the  resulting 
temperatures  for  each  box  are  compared  with  those  used  originally 
(block  6).  If  each  pair  of  temperatures  differ  by  no  more  than 
a specified  amount  (which  varies  with  latitude,  ranging  from  1.0K  near 
the  poles  to  0.2K  near  the  equator)  the  calculations  proceed  directly 
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to  block  7.  Otherwise,  the  initial  temperatures  are  adjusted  and 
the  sequence  of  blocks  2,  3,  and  5 repeated  (but  never  more  than 
six  times) . Notice  that  the  components  of  the  hydrologic  cycle 
are  not  recomputed  here . 

When  the  temperatures  have  converged  or  the  above  cycle 
has  been  completed  six  times,  the  temperatures  at  one  longitude 
(arbitrarily,  at  5°W)  are  compared  with  those  used  originally  to 
find  the  pressure  and  wind  fields  (block  7) . Whether  they  agree 
or  disagree  to  a specified  amount  (1.0K),  the  program  returns  to 
block  1.  With  non-convergence,  the  complete  cycle,  including  the 
calculation  of  the  components  of  the  hydrologic  cycle,  is  repeated, 
but  not  more  than  ten  times.  Satisfactory  convergence  is  usually 
achieved  after  no  more  than  four  cycles,  in  which  case  a final  pass 
is  made  through  block  5,  with  final  values  for  snow  and  ice  depths 
and  the  soil  moisture  content  made  along  the  way  (block  8) . After 
the  values  of  desired  variables  are  printed  out  (block  9)  the  program 
proceeds  to  the  next  month  (block  0) . At  the  end  of  each  year 
of  calculations  the  average  annual  temperature  in  each  box  and  the 
temperature  of  the  bottom  water  of  the  oceans  are  recomputed. 
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